Als nächstes wichtiges Mathe-Modul möchte ich hier an dieser Stelle noch eine schöne Routine für die p-q-Formel erarbeiten und einstellen.
Was steckt dahinter? Die Lösung einer quadratischen Gleichung in der Normalform x²+px+q=0 mit quadratischer Ergänzung nach x. Das Ergebnis ist die allgemeine pq Formel:
Und dafür bauen wir jetzt eine Routine, alle dafür benötigten Mittel haben wir bereits im Vorfeld erarbeitet.
... hat jetzt etwas länger gedauert, aber diese Assemblerroutine steht ja noch aus.
Die pq-Formel hat grundsätzlich 0,1 oder 2 Lösungen. Das muss man relativ am Anfang der Routine checken, sonst gibt es Chaos. Der erste Teil der Routine wird also das Ergebnis unter der Wurzel betrachten, nachdem es ausgerechnet worden ist: Wenn das Ergebnis negativ ist, gibt es insgesamt keine Lösung für x, wenn das Ergebnis 0 ist, dann gibt es genau eine Lösung für x und wenn das Ergebnis positiv ist, dann gibt es zwei Lösungen (x1 und x2). Also 1. Schritt das Ergebnis unter der Wurzel ausrechen und falls das Ergebnis größer 0 ist, die Wurzel davon ziehen.
Insgesamt nutzen wir für dieses Beispiel keine Fließkommazahlen, daher wird das Ergebnis immer nur eine Näherung sein. Allerdings brauchen wir von Anfang an auch negative Zahlen, daher beschäftigen wir uns zunächst damit: Da ich mit 8 Bit Zahlen rechnen möchte, nehme ich für die Kennung des Vorzeichens eine eine weitere Variable...
"Dumme" (oder nicht so "dumme") Frage: wofür wird die Formel benötigt? Was kann ich damit machen? Warum sollte ich mich damit beschäftigen?
https://mastodon.social/@cstrotm https://abbuc.social
Die pq-Formel ist wohl eine der essentiellen Formeln in der Mathematik. Damit löst man quadratische Gleichungen! Wenn man sich in der Mathematik bewegt, wird man wohl um die pq-Formel kaum herumkommen. Ob man die jetzt zwingend beim programmieren braucht, weiß ich auch nicht. Also zumindest beim Raycaster braucht es doch einiges an Mathematik (Pythagoras z.B.) - die pq-Formel ergänzt quasi diese kleine mathematische "Reise".....
... puh ... war doch aufwendiger, als gedacht! Aber noch in diesem Jahr fertig gestellt!! Hier nun ist der Code für die pq-Formel (wegen fehlender Fließkommazahlen natürlich nur eine Näherung):
Bitte gerne testen... hat mich gerade mein Hirn zum drehen gebracht 😉
LG und guten Rutsch Peter
var1 .by 08 ;p
var2 .by 01 ;q
var3 .by 00 ;x1
var4 .by 00 ;x2
minus1 .by 00 ;Minus für var1
minus2 .by 00 ;Minus für var2
minus3 .by 00 ;Minus für var3
minus4 .by 00 ;Minus für var4
temp1 .by 00,00 ;Zwischenspeicher1
temp2 .by 00,00 ;Zwischenspeicher2
anzahl .by 02 ;Anzahl der Lösungen
.proc main
;....................... pq-Formel .......
; x*x+px+q=0
; x1,2 = -p/2 ±√(p*p/4 - q)
; ----------------- p*p ausrechnen ----------------------------------------
lda var1 ; Lade p in den Akku
sta temp1 ; Speichere p im niedrigeren Byte von temp1
lda #0 ; Setze Akku auf 0 (wird als Akkumulator verwendet)
ldx #8 ; Initialisiere Schleifen-Zähler (8 Bits zu verarbeiten)
lsr temp1 ; Shift temp1 (niedriges Byte) um 1 nach rechts
lp1:
bcc lp2 ; Überspringe Addition, wenn kein Carry-Flag gesetzt
clc ; Clear Carry vor Addition
adc var1 ; Addiere var1 (aktuelles p) zum Akkumulator
lp2:
ror ; Rotate Right Akku (höheres Bit nach Carry)
ror temp1 ; Rotate Right temp1 (Carry ins höchstwertige Bit)
dex ; Dekrementiere Schleifen-Zähler
bne lp1 ; Wiederhole, solange Zähler nicht null
sta temp1+1 ; Speichere höheres Byte des Ergebnisses in temp1+1
; jetzt durch 4 teilen
clc
ror temp1+1
ror temp1
clc
ror temp1+1
ror temp1
;----------------- Subtraktion q -------------------------------------
#if .byte minus2=#0 ; für den Fall q ist positiv
sec
lda temp1
sbc var2
bcs weiter
#if .byte temp1+1=#0
lda #0
sta anzahl ;keine Lösung für x
jmp fertig
#end
dec temp1+1
jmp weiter
#end
clc ; für den Fall q ist negativ
lda temp1
adc var2
sta temp1
lda temp1+1
adc #0
sta temp1+1
lda temp1
weiter sta temp1
;-----------------------------------------
; Quadratwurzelberechnung (16-Bit)
; Eingabe: temp1/temp1+1 (16-Bit-Wert)
; Ausgabe: Quadratwurzel in temp2/temp2+1
;-----------------------------------------
lda #0 ; Akku auf 0 setzen
sta temp2 ; temp2 (niedriges Byte des Ergebnisses) initialisieren
sta temp2+1 ; temp2+1 (hohes Byte des Ergebnisses) initialisieren
ldx #8 ; Anzahl der Bits (8 Iterationen)
lp3
sec ; Carry-Flag setzen (für Subtraktion)
lda temp1+1 ; Höheres Byte von temp1 laden
sbc #64 ; Subtrahiere 64 (Schrittgröße)
tay ; Zwischenergebnis in Y speichern
lda temp2+1 ; Lade aktuelles Ergebnis (hohes Byte von temp2)
sbc temp2 ; Subtrahiere Zwischenergebnis
bcc lp4 ; Wenn Carry = 0, Rücknahme der Subtraktion
sty temp1+1 ; Speichere angepasstes Höheres Byte von temp1
sta temp2+1 ; Aktualisiere hohes Byte des Ergebnisses
lp4
rol temp2 ; Rotate temp2 (niedriges Byte) nach links
asl temp1 ; Shift temp1 (niedriges Byte) nach links
rol temp1+1 ; Shift temp1+1 (höheres Byte) nach links
rol temp2+1 ; Shift temp2+1 (hohes Byte des Ergebnisses) nach links
asl temp1 ; Wiederholen für Präzision
rol temp1+1
rol temp2+1
dex ; Dekrementiere Schleifen-Zähler
bne lp3 ; Wiederhole, bis alle Bits verarbeitet sind
;------------------- Wurzelberechnung fertig -----------------------------
lsr var1 ; p/2 bilden
#if .byte temp2=#0
lda #1
sta anzahl
sta minus3
mva var1 var3 ; p/2 als Ergebnis in x1 speichern
#if .byte minus1=#1
mva #0 minus3 ; in diesem Fall wird das Ergebnis positiv
#end
jmp fertig
#end
; ab hier jetzt x1 und x2 berechnen
; x1 = -p/2 + Wurzel / x2 = -p/2 - Wurzel
clc
lda var1 ; p/2
adc temp2 ; + Wurzelwert
sta var3 ; x1
lda var1 ; p/2
sec
sbc temp2 ; - Wurzelwert
bpl positiv ; den Wert als Betrag bilden
negativ
eor #$ff
clc
adc #1
positiv
sta var4 ; x2
; jetzt noch die Vorzeichen korrigieren
#if .byte temp2>var1
mva #1 minus3
mva #0 minus4
#if .byte minus1=#1 ; wenn p = minus
mva #0 minus3
mva #1 minus4
#end
jmp fertig
#end
lda #0
sta minus3
sta minus4
#if .byte minus1=#0
lda #1
sta minus3
sta minus4
#end
fertig
Erstmal danke für den Code, aber ich hätte da noch ein paar Anmerkungen: Sieh es bitte nicht als Kritik gegen Dich Peter, ich finde super gut, was Du hier etwas sendest.
Also:
sta temp1+1 ; Speichere höheres Byte des Ergebnisses in temp1+1
; jetzt durch 4 teilen
clc
ror temp1+1
ror temp1
clc
ror temp1+1
ror temp1
ist zu ersetzen durch:
lsr
ror temp1
lsr
ror temp1
sta temp1+1
weil kürzer und schneller 😉
IMHO sieht der Code scheußlich aus. Ich frage mich manchmal, wozu man Variablen überhaupt benennen kann. Warum benennst Du 'p' var1 und 'q' var2?
Zu was wird '#if .byte minus2=#0'. Das ist kein Assembler, auch weil der Mac65/Atasm das nicht übersetzt.
Welche Variablen darf/muss ich wie setzen?
q ist scheinbar auch nur 8bit breit.
Das Ergebnis landet wohl in x1,2. Wann genau gilt: Es gibt nichts, oder nur x1 oder x1 und x2. Das erschließt sich mir nicht so ganz.
Das scheint mir noch nicht ganz fertig zu sein.
JM2C
Einen hab ich noch: Dein Code enthält die Wurzelberechnung nach Lee Davison (http://6502org.wikidot.com/software-math-sqrt). Die Routinen von (https://www.codebase64.org/doku.php?id=base:fast_sqrt) sind etwas flotter. (so ca. 40%)
IMHO hatten wir das schon mal...
LG Lars
Hallo Lars,
ich wollte den Code unbedingt gestern noch reinstellen und vielleicht hätte es noch mehr Codetest und Bereinigung gebraucht. Um 01:15 konnte ich mich aber einfach nicht mehr so gut konzentrieren. Ich werde im neuen Jahr den Code noch überarbeiten, die #if / then ersetzen und die Variablen besser lesbar anpassen und noch mehr Erklärungen reinbringen. Ich schaue mir auch nochmal die andere Wurzelberechnung an... Danke für die Kritik. Der Code ist auch komplett frei und jeder darf ihn optimieren (wie du es schon beim teilen gemacht hast).
Liebe Grüße und jetzt aber erstmal mit gutem Rutsch ins neue Jahr kommen...
Peter
Hallo Ihr Lieben und ein frohes neues Jahr wünsche ich allen! Ich hoffe, alle sind gut ins neue Jahr hineingekommen.
Ich habe jetzt mal die Routine überarbeitet! Hoffe, daß es jetzt besser gefällt. Die langsamere Wurzeliteration habe ich noch so belassen, das will ich erst noch in meinen Code umsetzen und kann es dann hier einbauen.
p .by 00 ;p
q .by 00 ;q
x1 .by 00 ;x1
x2 .by 00 ;x2
minusp .by 00 ;Minus für p
minusq .by 00 ;Minus für q
minusx1 .by 00 ;Minus für x1
minusx2 .by 00 ;Minus für x2
temp1 .by 00,00 ;Zwischenspeicher1
temp2 .by 00,00 ;Zwischenspeicher2
anzahl .by 00 ;Anzahl der Lösungen
;....................... pq-Formel .......
; x*x+px+q=0
; x1,2 = -p/2 ±√(p*p/4 - q)
; p und q werden in diesem Beispiel bewußt auf 8 Bit
; beschränkt, da man ansonsten bei 16 Bit für p und q
; Zwischenspeicher und Wurzelfunktion auf 32 Bit hätte
; ausrechnen müssen - das würde hier den Rahmen sprengen
; ----------------- p*p ausrechnen ----------------------------------------
lda #8 ; 8-Bit-Wert in p laden
sta p
lda #0
sta minusp
lda #10 ; 8-Bit-Wert in q laden
sta q
lda #1
sta minusq
;--------------- Start pq-Routine -----------------------------------------
lda p ; Lade p in den Akku
sta temp1 ; Speichere p im niedrigeren Byte von temp1
lda #0 ; Setze Akku auf 0 (wird als Akkumulator verwendet)
ldx #8 ; Initialisiere Schleifen-Zähler (8 Bits zu verarbeiten)
lsr temp1 ; Shift temp1 (niedriges Byte) um 1 nach rechts
lp1
bcc lp2 ; Überspringe Addition, wenn kein Carry-Flag gesetzt
clc ; Clear Carry vor Addition
adc p ; Addiere p (aktuelles p) zum Akkumulator
lp2
ror ; Rotate Right Akku (höheres Bit nach Carry)
ror temp1 ; Rotate Right temp1 (Carry ins höchstwertige Bit)
dex ; Dekrementiere Schleifen-Zähler
bne lp1 ; Wiederhole, solange Zähler nicht null
; p*p durch 4 teilen
lsr
ror temp1
lsr
ror temp1
sta temp1+1
;----------------- Subtraktion q -------------------------------------
; für den Fall q ist positiv
lda minusq
cmp #0
bne sp2
sec
lda temp1
sbc q
bcs weiter
lda temp1+1
cmp #0
bne sp1
lda #0
sta anzahl ;keine Lösung für x
; hier Textausgabe: keine Lösung für x
jmp fertig
sp1
dec temp1+1
jmp weiter
sp2
; für den Fall q ist negativ
clc
lda temp1
adc q
sta temp1
lda temp1+1
adc #0
sta temp1+1
lda temp1
weiter sta temp1
;-----------------------------------------
; Quadratwurzelberechnung (16-Bit)
; Eingabe: temp1/temp1+1 (16-Bit-Wert)
; Ausgabe: Quadratwurzel in temp2/temp2+1
;-----------------------------------------
lda #0 ; Akku auf 0 setzen
sta temp2 ; temp2 (niedriges Byte des Ergebnisses) initialisieren
sta temp2+1 ; temp2+1 (hohes Byte des Ergebnisses) initialisieren
ldx #8 ; Anzahl der Bits (8 Iterationen)
lp3
sec ; Carry-Flag setzen (für Subtraktion)
lda temp1+1 ; Höheres Byte von temp1 laden
sbc #64 ; Subtrahiere 64 (Schrittgröße)
tay ; Zwischenergebnis in Y speichern
lda temp2+1 ; Lade aktuelles Ergebnis (hohes Byte von temp2)
sbc temp2 ; Subtrahiere Zwischenergebnis
bcc lp4 ; Wenn Carry = 0, Rücknahme der Subtraktion
sty temp1+1 ; Speichere angepasstes Höheres Byte von temp1
sta temp2+1 ; Aktualisiere hohes Byte des Ergebnisses
lp4
rol temp2 ; Rotate temp2 (niedriges Byte) nach links
asl temp1 ; Shift temp1 (niedriges Byte) nach links
rol temp1+1 ; Shift temp1+1 (höheres Byte) nach links
rol temp2+1 ; Shift temp2+1 (hohes Byte des Ergebnisses) nach links
asl temp1 ; Wiederholen für Präzision
rol temp1+1
rol temp2+1
dex ; Dekrementiere Schleifen-Zähler
bne lp3 ; Wiederhole, bis alle Bits verarbeitet sind
;------------------- Wurzelberechnung fertig -----------------------------
lsr p ; p/2 bilden
lda temp2
cmp #0
bne sp4
lda #1
sta anzahl
sta minusx1
lda p
sta x1 ; p/2 als Ergebnis in x1 speichern
lda minusp
cmp #1
bne sp3
lda #0
sta minusx1 ; in diesem Fall wird das Ergebnis positiv
sp3
; hier Textausgabe: "Lösung für x="x1
jmp fertig
sp4
; ab hier jetzt x1 und x2 berechnen
; x1 = -p/2 + Wurzel / x2 = -p/2 - Wurzel
clc
lda p ; p/2
adc temp2 ; + Wurzelwert
sta x1 ; x1
lda p ; p/2
sec
sbc temp2 ; - Wurzelwert
bpl positiv ; den Wert als Betrag bilden
negativ
eor #$ff
clc
adc #1
positiv
sta x2 ; x2
; jetzt noch die Vorzeichen korrigieren
lda p
cmp temp2
bcs sp6
lda #1
sta minusx1
lda #0
sta minusx2
lda minusp
cmp #1
bne sp5
lda #0
sta minusx1
lda #1
sta minusx2
sp5
jmp fertig
sp6
lda #0
sta minusx1
sta minusx2
lda minusp
cmp #0
bne sp7
lda #1
sta minusx1
sta minusx2
sp7
; hier Textausgabe: "Lösung für x1,2=" x1,x2
fertig
... anbei final die pq-Routine mit der schnellen Wurzelroutine:
p .by 00 ;p
q .by 00 ;q
x1 .by 00 ;x1
x2 .by 00 ;x2
minusp .by 00 ;Minus für p
minusq .by 00 ;Minus für q
minusx1 .by 00 ;Minus für x1
minusx2 .by 00 ;Minus für x2
temp .by 00,00
temp1 .by 00,00 ;Zwischenspeicher1
temp2 .by 00 ;Zwischenspeicher2
anzahl .by 00 ;Anzahl der Lösungen
;Iterationstabelle für die Wurzelberechnung:
iter: .byte $01,$02,$04,$08,$10,$20,$40,$80
;....................... pq-Formel .......
; x*x+px+q=0
; x1,2 = -p/2 ±√(p*p/4 - q)
; p und q werden in diesem Beispiel bewußt auf 8 Bit
; beschränkt, da man ansonsten bei 16 Bit für p und q
; Zwischenspeicher und Wurzelfunktion auf 32 Bit hätte
; ausrechnen müssen - das würde hier den Rahmen sprengen
; ----------------- p*p ausrechnen ----------------------------------------
lda #8 ; 8-Bit-Wert in p laden
sta p
lda #0
sta minusp
lda #10 ; 8-Bit-Wert in q laden
sta q
lda #1
sta minusq
;--------------- Start pq-Routine -----------------------------------------
lda p ; Lade p in den Akku
sta temp1 ; Speichere p im niedrigeren Byte von temp1
lda #0 ; Setze Akku auf 0 (wird als Akkumulator verwendet)
ldx #8 ; Initialisiere Schleifen-Zähler (8 Bits zu verarbeiten)
lsr temp1 ; Shift temp1 (niedriges Byte) um 1 nach rechts
lp1
bcc lp2 ; Überspringe Addition, wenn kein Carry-Flag gesetzt
clc ; Clear Carry vor Addition
adc p ; Addiere p (aktuelles p) zum Akkumulator
lp2
ror ; Rotate Right Akku (höheres Bit nach Carry)
ror temp1 ; Rotate Right temp1 (Carry ins höchstwertige Bit)
dex ; Dekrementiere Schleifen-Zähler
bne lp1 ; Wiederhole, solange Zähler nicht null
; p*p durch 4 teilen
lsr
ror temp1
lsr
ror temp1
sta temp1+1
;----------------- Subtraktion q -------------------------------------
; für den Fall q ist positiv
lda minusq
cmp #0
bne sp2
sec
lda temp1
sbc q
bcs weiter
lda temp1+1
cmp #0
bne sp1
lda #0
sta anzahl ;keine Lösung für x
; hier Textausgabe: keine Lösung für x
jmp fertig
sp1
dec temp1+1
jmp weiter
sp2
; für den Fall q ist negativ
clc
lda temp1
adc q
sta temp1
lda temp1+1
adc #0
sta temp1+1
lda temp1
weiter sta temp1
;-----------------------------------------
; Quadratwurzelberechnung (16-Bit)
; Eingabe: temp1/temp1+1 (16-Bit-Wert)
; Ausgabe: Quadratwurzel in temp2/temp2+1
;-----------------------------------------
ldy #$00
ldx #$07
clc
loopwu
tya
ora iter-1,X
sta temp+1
lda temp1+1
bcs wu0
cmp temp+1
bcc wu1 ; temp <= a
wu0
sbc temp+1
sta temp1+1
tya
ora iter,x
tay
wu1
asl temp1
rol temp1+1
dex
bne loopwu
; letzte Iteration
bcs wu2
sty temp+1
lda temp1
cmp #$80
lda temp1+1
sbc temp+1
bcc wu3
wu2
iny
wu3
sty temp2
;------------------- Wurzelberechnung fertig -----------------------------
lsr p ; p/2 bilden
lda temp2
cmp #0
bne sp4
lda #1
sta anzahl
sta minusx1
lda p
sta x1 ; p/2 als Ergebnis in x1 speichern
lda minusp
cmp #1
bne sp3
lda #0
sta minusx1 ; in diesem Fall wird das Ergebnis positiv
sp3
; hier Textausgabe: "Lösung für x="x1
jmp fertig
sp4
; ab hier jetzt x1 und x2 berechnen
; x1 = -p/2 + Wurzel / x2 = -p/2 - Wurzel
clc
lda p ; p/2
adc temp2 ; + Wurzelwert
sta x1 ; x1
lda p ; p/2
sec
sbc temp2 ; - Wurzelwert
bpl positiv ; den Wert als Betrag bilden
negativ
eor #$ff
clc
adc #1
positiv
sta x2 ; x2
; jetzt noch die Vorzeichen korrigieren
lda p
cmp temp2
bcs sp6
lda #1
sta minusx1
lda #0
sta minusx2
lda minusp
cmp #1
bne sp5
lda #0
sta minusx1
lda #1
sta minusx2
sp5
jmp fertig
sp6
lda #0
sta minusx1
sta minusx2
lda minusp
cmp #0
bne sp7
lda #1
sta minusx1
sta minusx2
sp7
; hier Textausgabe: "Lösung für x1,2=" x1,x2
fertig